seqblock2k <- function(object = "sbout2k.RData", id.vals, exact.vals = NULL, exact.restr = NULL, exact.alg = "single",	covar.vals = NULL, covar.restr = NULL, covars.ord = NULL, seed = NULL, assg.prob.stat = NULL, trim = NULL, assg.prob.method = NULL, assg.prob.kfac = NULL, assg.prob = NULL, file.name = NULL, query = FALSE, verbose = TRUE, ...){ 
	
	if(query == TRUE){
		object <- readline("Enter the name of the input file without quotation marks.  [E.g., sbout1.RData]  ")
	}
	
	load(object)  ## loads, doesn't export to wkspace

	## Rename object elements	
	x <- bdata$x  
	nid <- bdata$nid  ## names of id.vars
	nex <- bdata$nex  ## names of exact.vars
	ncv <- bdata$ncv  ## names of covar.vars
	rex <- bdata$rex  ## EXACT restricted values list
	rcv <- bdata$rcv  ## BLOCK restricted values list
	ocv <- bdata$ocv  ## block covariates order
	trn <- bdata$trn  ## treatment condition names
	apstat <- bdata$apstat ## assignment prob statistic
	mtrim <- bdata$mtrim ## trim fraction for apmeth=trimmean
	apmeth <- bdata$apmeth ## assignment prob method
	kfac <- bdata$kfac     ## assg multiple for method 'ktimes'
	assgpr <- bdata$assgpr ## assg probs for fixed prob
	orig <- bdata$orig  ## original unjittered data
	n.tr <- length(trn)
	
	## Define freq used quantities
	lid <- length(nid)
	lex <- length(nex)
	lcv <- length(ncv)	
	
	if(!(is.null(exact.restr))){
		rex <- exact.restr	
	}
	if(!(is.null(covar.restr))){
		rcv <- covar.restr	
	}
	if(!(is.null(covars.ord))){
		ocv <- covars.ord	
	}
	if(!(is.null(assg.prob.stat))){
		apstat <- assg.prob.stat
	}
	if(!(is.null(trim))){
		mtrim <- trim	
	}
	if(!(is.null(assg.prob.method))){
		apmeth <- assg.prob.method
	}	
	if(!(is.null(assg.prob.kfac))){
		kfac <- assg.prob.kfac
	}		
	if(!(is.null(assg.prob))){
		assgpr <- assg.prob
	}	

	if(query == TRUE){
		## ID
		id.vals <- NULL
		for(ii in 1:lid){
			tmp <- names(x)[ii]
			tmp2 <- readline(cat("Enter the value of '", tmp, "'.  ", sep="")) 
			id.vals <- append(id.vals, tmp2)
		}
		class(id.vals) <- class(x[,lid])  ## lid could be any of 1:lid
		
		if(sum(is.na(id.vals)) > 0){
			tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
			
			if(!(substr(tmp,1,1) %in% c("y", "n"))){
				ddd <- readline("The default is `yes'.  Continue? [y/n]  ")
					if(substr(ddd, 1, 1) != "n"){
						}else{
							tmp <- readline(cat("Warning: ID value(s)", which(is.na(id.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
						}
			}							
			if(substr(tmp,1,1) == "n"){
				stop()
			}
		}
		
		## Check for duplicated IDs if single id.var
		if((lid == 1) && (id.vals %in% x[, nid])){
			tmp <- readline(cat("Warning: ID value ", id.vals, " already used in the data.  Allow duplication and proceed? [y/n]", sep=""))
			if(substr(tmp,1,1) == "n"){
				stop()	
			}		
		}
		
		## EXACT
		exact.vals <- NULL
		if(lex > 0){
			for(ii in 1:lex){
        tmp <- nex[ii]
				tmp2 <- readline(cat("Enter the value of '", tmp, "'.  ", sep="")) 
				exact.vals <- append(exact.vals, tmp2)
			}
			class(exact.vals) <- class(x[,(lid+lex)]) ## could be any of (lid+1):(lid+lex)
		
			if(sum(is.na(exact.vals)) > 0){
				tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
				
				if(!(substr(tmp,1,1) %in% c("y", "n"))){
					ddd <- readline("The default is `yes'.  Continue? [y/n]  ")
					if(substr(ddd, 1, 1) != "n"){
						}else{
							tmp <- readline(cat("Warning: Exact blocking value(s)", which(is.na(exact.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
						}
					}							
				if(substr(tmp,1,1) == "n"){
					stop()
				}
			}
		}
		
		## EXACT ALGORITHM
		##   "single": if lex > 1, creates unique categories (4 for 2x2, e.g.)
		exact.alg <- "single"
		
		## COVARIATES
		covar.vals <- NULL
		if(lcv > 0){
			for(ii in 1:lcv){
				tmp <- ncv[ii]
				tmp2 <- readline(cat("Enter the value of '", tmp, "'.  ", sep="")) 
				covar.vals <- append(covar.vals, tmp2)
			}
			class(covar.vals) <- class(x[,lid+lex+lcv]) ## could be any of (lid+lex+1):(lid+lex+lcv)
			
			if(sum(is.na(covar.vals)) > 0){
				tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
				if(!(substr(tmp,1,1) %in% c("y", "n"))){
				ddd <- readline("The default is `yes'.  Continue? [y/n]  ")
					if(substr(ddd, 1, 1) != "n"){
						}else{
							tmp <- readline(cat("Warning: Blocking value(s)", which(is.na(covar.vals)), "is/are `NA'.  Proceed? [y/n] ", sep=" "))
						}
					}				
				if(substr(tmp,1,1) == "n"){
					stop()
				}
			}
		}
	
    ## SET RANDOM SEED:
		nnn <- readline("Would you like to specify a random seed for this experimental condition assignment?  [y/n]  ")
		
			if(!(substr(nnn,1,1) %in% c("y", "n"))){
				ddd <- readline("The default is `yes'.  Continue? [y/n]  ")
					if(substr(ddd, 1, 1) != "n"){
						}else{
							nnn <- readline(cat("Would you like to specify a random seed for this treatment assignment?  [y/n]  "))
						}
					}						
		if(substr(nnn,1,1) != "n"){
			seed <- as.numeric(readline("Enter the random seed.  "))
		}else{
			seed <- NULL
		}
		
		## OUTPUT FILE NAME:
		nnn <- readline("Would you like to specify the output data file name?  [y/n]  [If not, `sbout2k.RData' will be used.]  ")
    if(!(substr(nnn,1,1) %in% c("y", "n"))){
      ddd <- readline("The default is `yes'.  Continue? [y/n]  ")
      if(substr(ddd, 1, 1) != "n"){
      }else{
							nnn <- readline(cat("Would you like to specify the output data file name?  [y/n]  [If not, `sbout2k.RData' will be used.]  "))
      }
    }						
    if(substr(nnn,1,1) != "n"){
			file.name <- readline("Enter the file name.  Note: seqblock2k() expects a file of type `.RData'.  ")
    }else{
			file.name <- NULL
		}
  } ## close query==T 
	
	## Warn if duplicated ID values if single ID var
	if((lid==1) && (id.vals %in% x[, nid])){
			warning("Warning: ID value ", id.vals, " already used in the data.")
	}
	
	## check # exact vars vs. # exact vals
	if(lex != length(exact.vals)){
		stop(cat("Number of exact blocking variables (", lex, ") does not equal number of values (", length(exact.vals), ").\nSpecify exact block values and try again.\n", sep=""))
	} 
	
	## check if any exact val violates exact.restr
	if((lex > 0) && !(is.null(rex))){
		for(i in names(rex)){
			if(!(NA %in% rex[[i]])){
				rex[[i]] <- append(rex[[i]], NA)
			}
			wexrstr <- which(nex == i)
			if(!(exact.vals[wexrstr] %in% rex[[nex[wexrstr]]])){
				tmp <- readline(cat("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", i, "' is restricted.  Would you like to add ", exact.vals[wexrstr], " to the set?  [y/n]", sep=""))
				if(substr(tmp,1,1) != "n"){
					rex[[i]] <- append(rex[[i]], exact.vals[wexrstr])
				}else{
				stop("Exact blocking value number ", wexrstr, " (", exact.vals[wexrstr], ") is not in the set of values to which exact blocking variable `", nex[wexrstr], "' is restricted.  Respecify either the value, the restriction set, or both as needed." )
				}
			}
		}
	}
	
	## check if any block val violates covar.restr  
	if((lcv > 0) && (!is.null(rcv))){		for(i in names(rcv)){
			if(!(NA %in% rcv[[i]])){
				rcv[[i]] <- append(rcv[[i]], NA)
			}
			wcvrstr <- which(ncv == i)
			if(!(covar.vals[wcvrstr] %in% rcv[[ncv[wcvrstr]]])){
				tmp <- readline(cat("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", i, "' is restricted.  Would you like to add ", covar.vals[wcvrstr], " to the set?  [y/n]", sep=""))
				if(substr(tmp,1,1) != "n"){
					rcv[[i]] <- append(rcv[[i]], covar.vals[wcvrstr])
				}else{
				stop("Covariate blocking value number ", wcvrstr, " (", covar.vals[wcvrstr], ") is not in the set of values to which covariate blocking variable `", ncv[wcvrstr], "' is restricted.  Respecify either the value, the restriction set, or both as needed.")
				}
			}
		}
	}	
	
	if(!(length(ocv) == length(unique(ocv)))){
		stop("Prioritization of covariates for (number units considered) < (number covariates + 2) includes duplicate variable names.  Respecify 'covars.ord'.")	
	}	
	
	if(lex > 0){
		if(exact.alg == "single"){
			x.ex <- x[apply(t(t(x[,(lid+1):(lid+lex)])==exact.vals), 1, sum) == lex,]
		}
	}  
  
	## When no exact matches, use all previous data:
	if(is.null(exact.vals) || (nrow(x.ex)==0)){
		x.ex <- x
	} 
	
	prev <- nrow(x.ex)
	
	## Create data frame of unit 1 data
		## create id placeholders
	id.tmp <- rep(NA, lid)
	names(id.tmp) <- nid
		## create exact block placeholders
	ex.tmp <- rep(NA, lex)
	names(ex.tmp) <- nex

	qqq <- rbind(x.ex[,1:(ncol(x.ex)-1)], c(id.tmp, ex.tmp, covar.vals))
	qqq[nrow(qqq), 1:lid] <- id.vals
	if(!is.null(exact.vals)){
		qqq[nrow(qqq), (lid+1):(lid+lex)] <- exact.vals
	}
	
	orig <- rbind(orig, c(qqq[nrow(qqq),], Tr=NA))

	## Only when blocking covariates are used:
	## Default: covar order = order in dataframe
	if(lcv>0){
		if(is.null(ocv)){  
			ocv <- ncv  
		}
    
		## Select covars: if prev = 1, pick 1.  
		if(prev==1){
			cov.tmp <- ocv[1]  
   }
		## if prev>=2, pick prev-1 covars.
		##   to define MD, need num covars, lcv >= prev. 
		##   for distinct MD, need num covars, lcv >= prev+2
		if(prev >= 2){
			if(prev-2 >= lcv){
				cov.tmp <- ocv  ##changed to ocv from covars.ord
			}else{
        cov.tmp <- ocv[1:(prev-1)]
      } 
		}
		qqq.c <- data.frame(qqq[ , cov.tmp]) 
		names(qqq.c) <- cov.tmp
    
    ## If there is a variable with no variation, remove it from the vcov matrix for this assignment:
    wh.cut <- NULL
    for(jj in 1:(ncol(qqq.c))){
      if(isTRUE(var(qqq.c[, jj]) == 0)){
        wh.cut <- append(wh.cut, jj)
      }
    }
    if(length(wh.cut) > 0){
      n.qqq.c <- names(qqq.c)[-wh.cut]
      qqq.c <- data.frame(qqq.c[, -wh.cut])
      names(qqq.c) <- n.qqq.c
    }
    
    ## If new unit identical to old unit, jitter new unit's covar vals
		for(jj in 1:(nrow(qqq.c)-1)){
			if(isTRUE(sum((qqq.c[jj,] -qqq.c[nrow(qqq.c),])==0) ==length(qqq.c[1,]))){ 
				jit.val <- array()
				for(kk in 1:ncol(qqq.c)){
					jit.val[kk] <- jitter(qqq.c[nrow(qqq.c),kk])
				}
				qqq[nrow(qqq), cov.tmp] <- qqq.c[nrow(qqq.c),] <- jit.val 
			}
		}  ## To implement: jitter for identical rows w/ identical NAs.

	   ## replace NA's with variable means
	   for(jj in 1:ncol(qqq.c)){
	   	if(sum(is.na(qqq.c[ ,jj])) > 0 ){
	   		qqq.c[which(is.na(qqq.c[ ,jj])), jj] <- mean(qqq.c[ ,jj], na.rm=T)
	   	}
	   }

		mahmat <- mahal(qqq.c, cov(qqq.c))	
    
		tr.dist <- list()
		lrow <- mahmat[nrow(mahmat), 1:(ncol(mahmat)-1)]
	  ##	for(ii in levels(data1[, ncol(data1)])){}  ## good for factor TR
		for(ii in unique(x.ex[, ncol(x.ex)])){  ## good for character TR
			tr.dist[[ii]] <- lrow[which(x.ex[,ncol(x.ex)]==ii)]
		}      
	
		## Calculate summary for comparison
		if(apstat == "mean"){
			ms <- lapply(tr.dist, mean)  ## only prev assg'd tr's
		}
		if(apstat == "median"){
			ms <- lapply(tr.dist, median)  
		}	
		if(apstat == "trimmean"){
			ms <- lapply(tr.dist, mean, trim = mtrim)
		}
	
		## Sort previously assigned treatments, largest first
		tr.sort <- names(sort(unlist(ms), dec=T))  
		## Add treatments not yet assigned, in random order
		if(length(tr.sort) != n.tr){
			trn.unassg <- sample(trn[!(trn %in% names(ms))])
			nunassg <- length(trn.unassg)
			tr.sort <- c(trn.unassg, tr.sort)
		}

		if((length(ms) != length(tr.sort)) && (apmeth %in% c("prop", "prop2", "wprop"))){
			ms <- append(ms, rep(NA, nunassg), after=0)
			names(ms)[1:nunassg] <- trn.unassg
			## Give unassigned tr's double distance of max condition dist
			## (Quasi-minimization)
			ms[1:nunassg] <- 2*max(unlist(ms), na.rm=T)
		}
	
		if((length(tr.dist) != length(tr.sort)) && (apmeth == "wprop")){
			tr.dist <- append(tr.dist, rep(NA, nunassg), after=0)
			names(tr.dist)[1:nunassg] <- trn.unassg
			## Give unassigned tr's double distance of max condition dist, 
			##   with weight of one unit.
			## (Quasi-minimization)
			tr.dist[1:nunassg] <- max(unlist(ms), na.rm=T) 
		}
	}
	
	if(lcv == 0){ ## if there are only "exact" covariates used:
		tr.counts <- table(x.ex$Tr)
		trn.unassg <- trn[!(trn %in% c(x.ex$Tr))]
		if(length(trn.unassg) > 0){
			tr.counts <- c(tr.counts, rep(0, length(trn.unassg)))
			names(tr.counts)[(n.tr-length(trn.unassg)+1):n.tr] <- trn.unassg
		}
		p.lcv0 <- (1-tr.counts/sum(tr.counts))/sum(1-tr.counts/sum(tr.counts))
	}

	## Set assignment probability method
	if(apmeth == "ktimes"){
		p <- 	c(kfac/(kfac+n.tr-1), rep(1/(kfac+n.tr-1), n.tr-1))
	}
	
	if(apmeth=="fixed"){
		if(length(assgpr) != n.tr){
		stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) is length ", length(assgpr), ", while n.tr (number of treatment conditions) is ", n.tr, ".  Respecify and try again.")
		}		
		if(sum(assgpr) != 1){
		stop("assg.prob.method is `fixed', but assg.prob (vector of assignment probabilities) sums to ", sum(assgpr), " instead of 1.  Respecify and try again.")
		}
		p <- assgpr
	}		
	if(apmeth=="prop"){
		mssum <- sapply(ms, sum)
		p <- c(unname((mssum/sum(mssum))[tr.sort]))       
	}			 
	if(apmeth =="prop2"){  
		sumsq <- sapply(ms, sum)^2
		p <- c(unname((sumsq/sum(sumsq))[tr.sort]))
	}
	## Biases toward tr's w/fewer units
	if(apmeth=="wprop"){
		nnnj <- lapply(tr.dist, length)
		nnn <- sum(unlist(nnnj))
		div.func <- function(x){return(nnn/x)}
		nnnw <- lapply(nnnj, div.func)
		wdist <- unlist(ms)*unlist(nnnw)
		p <- c(unname((wdist/sum(wdist))[tr.sort]))
	}
	
	## set user-assigned seed
	if(!is.null(seed)){
		set.seed(seed)
	}

	## Assign new TR, given at least 1 "treated as continuous" covars
	if(lcv > 0){
    tr.new <- sample(tr.sort, 1, prob=p)
	}
  ## Assign new TR using correct order for probs and Tr condition names, EXACT covars only
	if(lcv == 0){
		tr.new <- sample(names(tr.counts), 1, prob=p.lcv0)
	}
  ## Assign new TR using same order for fixed probs and Tr condition names
	if(apmeth == "fixed"){
		tr.new <- sample(trn, 1, prob=p)
	}

	## Append to x
	x <- rbind(x, c(qqq[nrow(qqq),], Tr=tr.new))
	orig[nrow(orig),"Tr"] <- tr.new
	
	if(is.null(file.name)){
		file.name <- "sbout2k.RData"
	}
	
	if(substr(file.name, start=nchar(file.name)-5, stop=nchar(file.name)) !=".RData"){
		warning("Destination file name does not end in .RData.")
	}
	
	prev.dates <- bdata$datetime
	
	## Create storage list 'data2k'
	bdata <- list()
	bdata$x <- x
	bdata$nid <- nid
	bdata$nex <- nex
	bdata$ncv <- ncv
	bdata$rex <- rex  ## EXACT restricted values list
  if(!is.null(bdata$rex)){
    names(bdata$rex) <- nex    
  }
	bdata$rcv <- rcv  ## BLOCK restricted values list
  if(!is.null(bdata$rcv)){
    names(bdata$rcv) <- covar.vars
  }  
	bdata$ocv <- ocv  ## block covariates order
	bdata$trn <- trn
	bdata$apstat <- apstat ## assignment prob statistic
	bdata$mtrim <- mtrim ## trim fraction for apmeth=trimmean
	bdata$apmeth <- apmeth ## assignment prob method
	bdata$kfac <- kfac     ## assg probs for method=ktimes
	bdata$assgpr <- assgpr ## assg probs for fixed prob
	if(lcv>0){
		bdata$trd <- tr.dist
    bdata$tr.sort <- tr.sort
    bdata$p <- p
	}else{
		bdata$trcount <- tr.counts
	}
	bdata$datetime <- append(prev.dates, date())
	bdata$orig <- orig
	
	save(bdata, file = file.name)
	if(verbose == TRUE){
		cat("Unit 1:", nrow(x), " data stored as file ", file.name, ".\nThe current working directory is ", getwd(), "\n", sep="")
		cat("Unit ", nrow(x), " assigned to ", tr.new, ".\n", sep="")
		cat("The new data as entered:\n")
		print(x[nrow(x),])
	}	
	return(bdata)
}
